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Abstract. We compute the warping of a disc induced by an inclined dipole. We consider a magnetised star surrounded 
by a thin Keplerian diamagnetic disc with an inner edge that corotates with the star. We suppose the stellar field is 
a dipole with an axis that is slightly misaligned with the stellar rotation axis. The rotation axes of the disc material 
orbiting at large distances from the star and that of the star are supposed to coincide. The misalignment of the 
magnetic and rotation axes results in the magnetic pressure not being the same on the upper and lower surfaces of the 
disc. The resultant net vertical force produces a warp which appears stationary in a frame corotating with the star. 
We find that, if viscosity is large enough (a ~ 0.01-0.1) to damp bending waves as they propagate away, a smoothly 
varying warp of the inner region of the disc is produced. The amplitude of the warp can easily be on the order of ten 
percent of the disc inner radius for reasonably small misalignment angles (less than 30 degrees). Viscous damping also 
introduces a phase shift between the warp and the forcing torque, which results in the locations of maximum elevation 
above the disc forming a trailing spiral pattern. We apply these results to recent observations of AA Tau, and show 
that the variability of its light curve, which occurs with a period comparable to the expected stellar rotation period, 
could be due to obscuration produced by a warp configuration of the type we obtain. 
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1. Introduction 

Objects accreting material through an accretion disc very commonly contain a significant magnetic field. This is the 
case for accreting white dwarfs in cataclysmic variables, some X-ray binary pulsars and at least some classical T Tauri 
stars (CCTS). 



It was first suggested by Bertout et al. ( 1988 ), as a result of the detection of bright stellar spots, that CTTS may 
accrete along stellar magnetic field lines. This picture has been further supported by a wide array of observational 
evidence (see Najita et al. |2000 an d refe rences therein), i nclud ing spectroscopic indications of infalling materia l onto 



the stellar surface (Edwards et al. 1994; Hartmann et al. 1994) and the low spin rate of CCTS (Bouvier et al. 1993; 
Edwards et al. |1993| ). 

Since T Tauri stars have a large convective envelope, it is likely that at least part of their magnetic field is generated 
through a dynamo process. However, there may also be a fossil component originating from the molecular cloud out 
of which the star formed (Tayler 1987). Recent Zeeman measurements indicate relatively strong field strength at the 
surface of T Tauri stars, on the order of one kilogauss (Guenther et al. 199£; Johns-Krull et al. 



1999) 



It is not known 

what the structure of the field is. At some distance from the star the dipolar component probably dominates, but 
whether this is the case in the magnetosphere is not clear. However, observations cannot rule out such a coherent 



field structure (Montmerle et al. 1994), and numerical simulations of nonlinear stellar dynamos indicate that a steady 
dipole mode is the most easily excited one (Brandenburg et al. 1989] ). 
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Interaction between the stellar magnetic field and the accretion disc has very important consequences for the disc 
structure, the accretion process (see Ghosh & Lamb 1991 and references therein) and the evolution of the stellar 
rotation (Konigl 1991). In particular, the disc is truncated by the magnetic torque, so that it does not extend down 
to the stellar surface (Gosh & Lamb 197£). The location of the disc inner radius is determined by the condition that 
magnetic and viscous torques balance. For CTTS, the radius of the inner cavity is believed to be a few stellar radii 
(see, e.g., Wang 1995| ) 



So far, there are only a few numerical simulations of disc-stellar magnetic field interactions (Hayashi et al. 1996; 
Miller & Stone 1997; Goodson et al. 1997; Kudoh et al. 199E). They all show the disc-magnetosphere interaction 



to be complex and sensitive to initial and boundary conditions. At this stage, it is not clear what final form a full 
theoretical model is likely to take. However, analytical or semi-analytical simplified models can still be valuable in 
pointing out some important processes that may arise in these systems, and the goal of this paper is to describe one 
of these processes. 

We note that the magnetic axis and the rotation axis of the disc at large distances from the star may not be 
aligned, although often, for simplicity, they are assumed to coincide. We assume here that the stellar rotation axis 
and the disc rotation axis at large distances coincide. Misalignment would then occur if, for instance, the star were to 
generate a dipole field with magnetic axis misaligned with its spin axis (like in the case of the Earth) . It is possible 
that the interaction with the disc would lead to some evolution of the misalignment angle, but the details are likely 
to depend on the processes which generate the field. In any case, when such a misalignment is present, the magnetic 
pressure is not the same on the upper and lower surfaces of the disc. This mismatch generates a net vertical force 
which excites bending waves and warps inner parts of the disc (Aly |1980 ). 

Bending instabilities in a disc subject to a stellar magnetic dipole have been investigated by Agapitou et al. ( |1997| 
hereafter APT). APT calculated the global bending modes of a disc permeated by both an internally produced poloidal 
magnetic field and an external dipole field with axis aligned with the disc rotation axis (in this case no warp is induced 
by the dipole configuration, but free bending modes can be excited by a perturbation which takes the disc out of its 
equilibrium plane). They found that instability could occur if the magnetic and centrifugal forces were comparable in 
some region of the disc. They pointed out that such instabilities may result in the periodic variability observed in the 
light curve of many CTTS. 

Lai Q1999D studied the warping of a disc induced by an inclined dipole. He calculated the magnetic torque exerted 
by an inclined dipole on a disc, and studied the stability against vertical displacements of a disc subject to such a 
torque. In the terms of the APT analysis, he studied the stability of low frequency (as measured in an inertial frame) 
bending modes corresponding to the modified tilt mode as discussed in APT. We note that, when considering the 
structure of the disc subject to the inclined dipole, he did not take into account the effects of the distortion of the 
disc itself on its response, which can have important consequences on the dynamics through wave propagation. But 
he added the effects of a toroidal field, assumed to be generated by winding up a penetrating vertical field, on the 
magnetic pressure determining the vertical force on the disc. This contribution is phase shifted with respect to the 
other contributions and may thus (if not counteracted) cause the modified tilt mode to become unstable, resulting in 
spontaneous warping. To decide whether this mode can be destabilised requires detailed consideration of the effects 
of wave propagation and viscosity. We comment that un der so me conditions warps diffuse away on a timescale much 
shorter than the viscous timesc ale (P apaloizou & Pringle 1983 ) or propagate away with a velocity on the order of the 
sound speed (Papaloizou & Lin 1995 ) resulting in stabilisation. 

In this paper, we calculate the structure of a thin Keplcrian disc subject to an inclined dipole, taking into account 
the effects of the distortion of the disc itself on its response, i.e. the full dynamics of the system. For simplicity, 
we suppose that the disc is diamagnetic, so that it is not permeated by the external stellar field. In principle, the 
calculations presented here could be extended to more general cases. However, if the disc were not diamagnetic, 
wrapping of field lines would probably become important, leading to the possible disruption of the magnetosphcre 
(see, e.g., Mikic & Linker 1994 ). Also we do not address here the physical processes of accretion or plasma entry into 
the stellar magnetosphere. We note that because the warp induced in the inner disc appears steady in a frame rotating 
with the star, any resulting variability would have the same period as that of the star. 

We comment that the generation of spontaneous warping does not apply to the calculations we present here, since 
we study a response which is forced by the inclined dipo le and has a pattern speed equal to the rotation rate of the 
star. Thus, in contrast to the considerations of Lai ( |l999| ), it is not a modified tilt mode. 

This work has been motivated by a recent study of Bouvier et al. (199E) who report that the light curve of the 
CTTS A A Tau displays photometric, spectroscopic and polarimetric variations on timescales from a few hours to 
several weeks. The most striking feature of this light curve is a photometric variability with a period comparable 
to the expected rotation period of the star. This has been interpreted by Bouvier et al. ( 1999 ) as being due to the 
occultation of the star by a warp of the inner disc (the system is observed almost edge-on) . The authors speculated 



Terquem & Papaloizou: The response of an accretion disc to an inclined dipole 



3 



that the warp could be produced by an inclined dipole. We note that Bouvier et al. (1999) did not consider AA Tau 
as being a special case as far as its properties are concerned. They pointed out that only its light curve is unusual, 
and they interpreted it as being due to the fact that the system is seen almost edge-on. In other words, warping of 
the inner parts of CTTS discs would not be uncommon, but it could be seen only for particular viewing angles. 

The plan of the paper is as follows: We begin by considering an equilibrium configuration where the axis of the 
dipole and the rotation axis of the disc are aligned. This is described in § |^. We then perturb this equilibrium by 
slightly inclining the dipole. In § [s]we calculate the resulting perturbed magnetic field and derive the integro-differential 
equation which has to be solved for the disc vertical displacement. This equation is solved in the WKB approximation, 
which is valid when the wavelength of the bending waves excited in response to the perturbation is small compared to 
the disc radius. To allow these waves to damp as they propagate away (and therefore to damp very small wavelength 
oscillations, which would be unphysical), we include the effects of viscous damping in the integro-differential equation. 
In § ^ we so lve this e quatio n numerically and present the results for two different magnetic field equilibria (derived 



by Aly 1980 and Low 1986 ). In both cases we find that, if the viscosity is large enough to damp the waves as they 
propagate away, a smooth warp configuration of the disc inner parts can exist. The elevation above the disc equilibrium 
plane (i.e. the stellar equatorial plane) can easily be on the order of ten percent of the disc inner radius for reasonably 
small misalignment angles (less than 30 degrees). If the viscosity is too small to damp the waves efficiently, the disc 
inner parts may be disrupted, and it is likely that the evolution is then highly time-dependent. In § ^ we apply these 
results to the case of AA Tau, and show that the variability of its light curve, which occurs with a period expected to 
be the stellar rotation period, can plausibly be explained by a warp configuration of the type we obtain. 

2. Disc and aligned dipole 

We begin by considering a thin disc configuration such that the gas orbits a central rotating star with a dipole field 
and where the magnetic axis and rotation axes are all aligned. In this situation, the configuration is axisymmetric. 
For convenience, we work in a frame corotating with the central star with angular velocity oj. The disc is truncated in 
its inner parts by the magnetic torque and we will suppose that the inner edge corotates with the star. 



2.1. Magnetic field 

The dipole field B ea;t due to the central star induces azimuthal currents in the conducting disc which in turn generate 
an additional poloidal magnetic field B f j. 

Here, as in APT, we shall assume that the field external to the disc and central star can be approximated as curl 
free, i.e. assumed to be a vacuum field (the currents external to the disc and star are at infinity). 

The total axisymmetric poloidal magnetic field, B e:E f + B^, is noted B = (B r , 0, B z ), where we use cylindrical polar 
coordinates (r, cp, z). The associated Cartesian coordinates are (x, y, z). The flux function ip is such that: 

-1 dip 1 dip 

B r = — — a,n&B z = - — . (1) 
r oz r or 

For the stellar dipole field: 

/ _ Var 1 

Vext ~ ( r 2 + 2)3/2' y z > 

where the magnitude of the stellar magnetic dipole moment is fid = B*R^, with B* and i?* being the stellar magnetic 
field and radius, respectively. 

For the general axisymmetric poloidal field, the associated current density is j = (0,j v ,0). For an infinitcsimally thin 
disc, as we consider here, we define the vertically integrated azimuthal component of the current density, J, such that: 



J = j v dz. 

J — OO 

By integrating the azimuthal component of Ampere's law through the disc, we obtain: 

B+ = J/2, (3) 

where Bf denotes the radial component of the magnetic field just outside the upper surface of the disc. Throughout 
this paper we use MKSA units, and /io denotes the permeability of the vacuum. B r is antisymmetric with respect to 
reflection in the disc mid-plane so that its value just outside the lower surface of the disc is B~ = —Bf. 
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2.2. Force balance 

Neglecting pressure forces as in APT, the vertical integration of the radial component of the momentum equation 
yields the condition for radial equilibrium as: 

s|5 = Erf} 2 + JB Z , (4) 

where £ is the surface density, fl is the disc angular velocity measured in an inertial frame, and $ is the gravitational 
potential here taken to be due to a central point mass, M*, such that: 



Vr 2 + z 2 



If B z 7^ 0, in some regions of the disc a variety of configurations are possible (see APT). These include cases where 
inner field lines cross the disc and join to the central star and thus may be assumed to corotate with it. Outer field 
lines may be open in the case of an infinite disc or a finite disc with appropriate boundary conditions. 

In this paper, for simplicity, we shall perform calculations for the special case where the field is excluded from the 
disc. This situation arises when the disc is perfectly conducting. Then only surface currents flow in the disc and they 
screen the stellar dipole field from the disc interior, i.e. they produce an additional field which cancels TS ext in the 
disc interior. We note that, since the vertical component of the field is continuous at the disc surface, B z just outside 
the disc is zero in this case. The same does not apply for B r , and in general and B~ are non zero. 

We remark that, since there is no Lorentz force acting in the disc when B z = 0, the angular velocity given by 
equation (Q) is Keplerian. 

3. Disc and slightly misaligned dipole 

We suppose that the system is perturbed from the axisymmetric equilibrium state described above by introducing a 
slight misalignment between the magnetic axis of the central dipole and the rotation axis of the disc and star (the z 
axis). This produces a non axisymmetric response in the disc that can be described using linear perturbation theory. 



The main features of this response is that it takes the form of a warping of the disc, as indicated by Aly (1980), 
together with the additional feature of the excitation of bending waves. The response is naturally largest in the inner 
parts of the disc where the magnetic field is strongest. As we indicate by considering specific examples in § |], it can 
take the form of a steeply changing inclination of the inner disc orbits which can make it appear to have an inner wall. 

To calculate the geometry of the disc (that is its elevation above its equilibrium plane) , we extend the analysis of 
APT. To compute the free bending modes of the disc, APT solved an eigenvalue problem where the eigenvalues were 
the mode frequencies and the eigenfunctions their amplitudes. Here we are interested in a response with frequency 
equal to that of the forcing term. Since the dipole is anchored in the star, this is the angular velocity of the star. 
The amplitude of the mode (its spatial dependence) can be found from the mode equation of APT with a specified 
frequency and the addition of the forcing term. However, a different equilibrium field is adopted, since here, in contrast 
to APT, there is no internally produced field in the disc. 

3.1. Perturbed magnetic field 

To calculate the response we suppose the central dipole moment is rotated in the (x, z) plane through a small angle 8 
being the inclination to the z axis. The dipole moment is then given by: 

= (Md£,0=Md)- ( 5 ) 
This contributes to a potential: 

*ke*t = ^ ( 6 ) 

where r denotes the position vector, which produces a radial external magnetic field perturbation B' r ext = d<&' M ext /dr. 
Just outside the disc surfaces, this is given by the real part of: 

R /+ _ R /- _ g/frfg expjitp) 

n r,ext ~ n r.ext ~ „q \' ) 
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Thus the problem reduces to the calculation of the disc response to a field perturbation with azimuthal mode 
number m = 1. This, when acting with the unperturbed azimuthal current, produces a vertical Lorentz force which 
tends to warp the disc. In other terms, when the dipole is misaligned, the magnetic pressure force is not the same on 
the upper and lower surfaces of the disc. This mismatch generates a net vertical pressure force which tends to warp 
the disc. 

We thus introduce the Lagrangian displacement £ which, in a razor-thin disc approximation, has the form: 

€ = (0,0,6,)- 

The only non- negligible component is the vertical one (APT), and £ z represents the elevation above the disc equilibrium 
plane. In keeping with the form of the external magnetic field perturbation, we make all the perturbed quantities 
complex by taking their (/j-dependence to be through a factor exp(i^), which henceforth will be dropped. The physical 
perturbations will be recovered by taking the real part of these complex quantities. The Eulerian perturbations of the 
various quantities are denoted by a prime. 

The perturbation of the magnetic field interior to the disc, B', is related to 6 by the flux freezing condition: 

B' = (5;,5;,B^ = VX(«XB). (8) 

The non-zero components of B' take the form: 

B r = , and B z = . (9) 

dz r or 

Since B r is antisymmetric with respect to reflection in the disc mid-plane and £ 2 is independent of z in the thin disc 
approximation, B' z is antisymmetric and B' r is symmetric. 

As in APT, the vertical component of the perturbed field must be matched to the vertical component of a perturbed 
vacuum field exterior to the disc. Here this is taken to be a potential field. Therefore we have: 

^=B' Z \ (10) 

where B'+ is the value of the vertical field perturbation just outside the disc surface and Q' M is the magnetic potential 
associated with the external field perturbation. To find §' M , we first subtract out the contribution of the external field 
perturbation arising from the tilted dipole. At the disc surface, this is (eq. ||): 

$' fill 

where as above the factor exp(icp) has been dropped. With this contribution removed, the residual potential has no 
singularity outside the disc. It can be calculated from equation ( |To| ) in an analogous manner to finding the gravitational 
potential due to a disc surface density distribution, with 27rGS (where G is the gravitational constant) being replaced 
by B' z + (see Tagger et al. 199C; Spruit et al. 1995; APT). Following this procedure, $' M may be written as the sum of 



a Poisson integral and &' M ext : 

$M = $M,d + $ M,e*t, (12) 

where 

& M d = -— r° r B ^^° os ^ r ' dr ' d ^ ( i3) 

27r J R . J Q y/ r '2 + r 2 ~ 2rr' cos(tp') + z 2 ' 

where Ri and R a are the inner and outer radii of the disc, respectively. In the above integral, the (/3-dependence of the 
perturbed field has been taken into account. Here again, the factor exp(iy), to which <f>' M d is proportional, has been 
dropped. 

Since the disc is perfectly conducting, the vertical component of the field given by equation (^) is continuous at 
the disc surface. Therefore: 

r dr [ ' 

This expression for B' z + can be used in the above integral. 
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The radial component of the magnetic held just outside the surfaces of the disc is given by: 



where equations ( |lT| ) and ([L2|) have been used to write the last expression (we have approximated the derivatives at 
the disc surfaces by their value at z — because the disc is inhnitesimally thin). 

3.2. The disc vertical displacement 

The vertical component of the perturbed Lorentz force integrated vertically through the disc is (see APT): 
2R+ R'+ 

F > dz = _^r£^ (16) 
) MO 

The term prop orti onal to dB z /dr, which was present in APT, is zero here since B z = in the disc interior and at its 



surfaces (see § |2.2| ) 



We note that this force is simply the perturbed magnetic pressure force P' m vertically integrated through the disc. 
The pressure force vertically integrated through the disc is indeed: 

m ~ 2 M o + 2 Mo ' (17j 
(B z = at the disc surfaces), so that: 

pL = _BlB* + B £ B^ = _2_BlB£ > (18) 
Mo Mo Mo 

where we have used the symmetries of the field. We note that P' m is non zero because the perturbation induced by 
the misaligned dipole has an opposite effect on the upper and lower surfaces of the disc: it increases the magnitude of 
the radial component of the magnetic field on one of the surfaces from \B^\ to \B+\ + \B'+\ while decreasing it on the 
other surface from to \B+\ - \B'+\. 

The vertically integrated z-component of the perturbed equation of motion is: 

= 6+/ Ffr (19) 



Di 2 



2 = 



where D/Di denotes the convective derivative. 

In the problem considered here, the solution is steady in a frame rotating with the central star angular velocity w. 
Then: 

£^ = -£(c-fi) 2 £ 2 . (20) 

Also, for a point mass potential: 
<9 2 $\ GM 



., . , - ; =K, (21) 

where J7k is the Keplerian angular velocity. Using equations (|l5|), (p^|), ( pc| ) and (|l|), equation (|l^) becomes: 

2=0 r3 



Or 



(22) 



The potential <fr' M d can be expressed in term of £ z and its derivative with respect to r using cq. (|13|) and (|lj) . Therefore 
equation ( |22| ) gives a linear equation for determining the vertical displacement £ 2 induced by the tilted external dipole. 
The term proportional to [id in this equation acts like a forcing term for the vertical displacement and it causes the 
disc to become warped and may also excite bending waves. 



As noted above, since there is no Lorentz force acting in the disc interior, f2 in equation (22) is the Keplerian angular 
velocity. 
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We remark that equation (|22|) is the same as equation (14) of APT with B z — 0, but with an additional nonhomo- 
geneous term proportional to /i^. As noted above, this is expected as here we calculate the response which is forced by 
the nonaligned dipole. We remark that this nonhomogeneous term is simply the pressure force due to the misaligned 
component of the dipole vertically integrated through the disc: 



P' 



2B+B'+ ext _ 2B+ 2^6 



Mo Mo 

where we have used equation (Q) (see also eq. [24] of Aly 1980 ). 



(23) 



3.3. WKB waves 

When [id is set to zero in equation (p2|), this equation has solutions corresponding to free bending waves. In the local 
limit, these can be found by assuming that £ 2 cx exp(ifcr), where k is the radial wavenumber, assumed to be very large 
compared to 1/r. The integral in equation ( |l3| ) can be evaluated in the WKB approximation to give, after using (14) 
(see APT and references therein): 



M 



\k\ 



-ikB+£ z 
\k\ 



The local dispersion relation derived from equation (B2J) is then 

(cj-n K f = n 2 K + ^-\k\, 



(24) 



(25) 



where we have used the fact that the angular velocity is Keplerian. We see from equation ( j25| ) that bending waves 
with m — I propagate (with \k\ > 0) exterior to the Lindblad resonances where (w — Q^) — ±£!k- In this paper, we 
are concerned with the situation where the disc is terminated at an inner boundary (r = Ri) where the local Keplerian 
rotation rate is close to corotation with the central star, i.e. flK{Ri) = w. In such a case, only the outer Lindblad 
resonance (OLR), where J7k = will exist within the disc, and it occurs at a radius r — 1.59-Rj. 
Significantly beyond the OLR, the wavenumber associated with the waves is given by: 



^ w 2 r£ r 3 
1 1 2(B+) 2 13 1 



(26) 



where we define j3 such that it would be the ratio of magnetic to centrifugal forces if we had a vertical field B z ~ £>+ 
(see § Although B z — here, we use f3 as a measure of the strength of This is expected to be at most of order 
unity at the inner edge of the disc and then to decrease rapidly outwards, as the magnetic field decreases (see § |J). 

If (3 decreases as r increases, the forcing term proportional to \id in equation ( ^2|) will excite bending waves in the 
disc that propagate outwards with increasing wavenumber until they are damped, much as in the case of Saturn's 
rings (see Goldreich & Tremaine 1978 ; Shu 1984 ). In some cases, the wavelength may be so short at the OLR that 
dissipative processes prevent waves from being properly launched (see below). 



3.4- Viscous dissipation 

In order to allow the waves to damp as they propagate away from the OLR, we add an additional viscous force per unit 
mass to the right hand side of equation (^) of the form i>d 2 v' z / dr 2 , where v' z is the Eulerian perturbed vertical velocity 
and v is a kinematic viscosity. Neglecting the variation of Q k, which is reasonable for short wavelength disturbances, 
this can also be written as w(51k — uj)d 2 £, z /dr 2 . For v we use a standar d 'a' prescription such that v — aH 2 ilK, where 
H is the disc semi-thickness and a is a constant (Shakura & Sunyaev 1973 ). 

With the incorporation of the above viscous force, equation (E2h, which gives the forced response of the disc, becomes: 



-(ui - ft) 2 + n|] £ 2 - iv(fl K - lo 



dr 2 



2B+ 
MoS 



9$ 



M.d 



dr 



2fi d 5 



(27) 



We now comment briefly on the possible origin of this a viscosity. So far, the only process which has been shown 
to initiate and sustain turbulence in Keplerian discs is the magnetohydrodynamic Balbus-Hawley instability, and it 
does lend itself to an a formalism (Balbus & Hawley 1998 and references therein). For the disc we consider here, which 
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is not permeated by the external magnetic field, it is difficult to justify the existence of an a viscosity. However, in 
reality the disc is probably permeated to at least a small extent by the external field (Ghosh & Lamb 1979). If the 
ratio of magnetic to thermal pressure is smaller than unity, the disc is then subject to the BH instability, the internal 
disc field is amplified and bending waves are damped. The values of a produced by the BH instability range roughly 
from 10~ 3 to 0.1 (Hawley et al. 1995; Brandenburg et al. 1995), the largest values corresponding to the case where 



the magnetic field varies on a scale large compared to the disc semithickness (Hawley |2000 ) 



4. Numerical response calculations 

4--1- Unperturbed magnetic field 

For the model calculations presented here, we adopt a radial magnetic field corresponding to the situation when the 



centr al dip ole flux is completely excluded from the disc. We consider both the solutions computed by Aly ( 198C ) and 
Low (|1986|) . Aly's solution corresponds to the case when all the dipole flux goes through the magnetospheric cavity in 
the middle of the disc, and is singular at the disc inner edge. Low gives a solution that is non singular with some flux 
escaping to infinity. 

In these cases we have, for r larger than some radius r\e: 



s-r = Tj Bo( — 



r 



2 

^) - I 



-| ±1/2 

(28) 



where we have defined: 



* = %■ < 29 » 

In equation (|28|), the upper sign corresponds to Low's solution, whereas the lower sign corresponds to Aly's solution. 
Of course, in both cases B z = since the disc is diamagnetic. Note that the radial field at the disc surface has opposite 
sign in the two cases. 

We comment that -B+ is negative and positive for Low's and Aly's cases, respectively, which is the opposite of 
the solution obtained by these authors. This is because we have arbitrarily chosen B z to be positive when the dipole 
moment fi d is positive, in contrast to these authors. 

In a two dimensional model, rg is the disc inner radius. However, more realistically it is expected to differ from 
this by an amount comparable to the vertical thickness of the disc. We adopted the procedure of taking the disc inner 
radius Ri to be slightly larger than rs, and to check convergence of the results by decreasing (Ri — re). For Aly's 
solution, convergence requires that the disc surface density increases rapidly enough towards Ri such as to limit the 
Lorentz force there. 



4-2. Results 

We have performed global normal mode calculations for disc models with the unperturbed magnetic field described 
above. 

We solve equation ( |27| ) considered as an integro-differential equation for the response function £ z by the method 
described in APT. We use the dimensionless radius 137 = r/R a , so that the outer radius is va = 1. In these units, the 
disc inner radius on the computational grid is taken to be Wi = Ri/R =0.1, and to corotate with the central star so 
that oj = (Ri)- 

The radial interval [wi, 1] is divided into a grid of n r equally spaced points at positions (wi) I=1 with a spacing 
Aw 1 = vji+i — w 1. Equation ( |27j ) is solved as in APT by discretizing it on the grid so converting it into a matrix 
inversion problem for (£ z ) /=1 . In the calculations presented below, n r was varied between 500 and 700, but 300 
already gave satisfactory results. 



As already noted in § 3.2, we use the magnetic support /3, which would be the ratio of the Lorentz force to the 
centrifugal force in the disc if B z were non zero and on the order of _B+ : 

_ 2(g+)V 
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where G is the gravitational constant. Although B z = in the disc, (3 will be used as a measure of the strength of the 
field . Wc found it convenient to define the dimcnsionlcss quantities E = E/Ej, where Ej if the surface density at 
the disc inner radius, B+ = B+ /B , and (3 — (3 /(3a with: 

O = -tgkr- (31) 

We then have: 

- 2H7 2 (B+f , 
0= y^-- (32) 

For both Low's and Aly's magnetic field, we consider models where the magnetic support is large (reaching values 
of order unity) close to the disc inner edge and decreases rapidly with radius. For Aly's solution, we take f3 = 1 at the 
disc inner edge. For Low's solution, j3 has to vanish at the same location as B^ , i.e. at wb = Tb/Ro, otherwise £ z is 
infinite there. We then take /3 = at wb and /3 = 1 at the location where B+ has its maximum. From the value of (3 
and _B+ at the disc inner radius, we can calculate [3$. 

From the profile of f3, we calculate E using equation (|32"|). We actually take the profile of E which gives the desired 
f3 near the disc inner edge, and we then match it on to a r~ 3 / 2 profile further away in the disc (since B^ decreases 
rapidly with radius, we do not need to worry about the profile of E for w larger than a value which turns out to be 
about 0.2). 

We note that our results do not depend on the numerical values of E^, Bq, and Ri (or equivalently R ) taken 
individually, but only on (3q and 5. For Ri we could take a few times and then compute R Q = lORi. We could 
then get a value of -Bo by using equation (p9| ) with fid = B r R^. Then E^ would follow from the expression (|3|) for [3q. 
Depending on the value of Ej , the amount of mass in the annulus between Ri and R a could be modified by changing 
the profile of E for w > 0.2. 

An external poloidal magnetic field tends to squeeze the disc. Therefore, we choose a form of the aspect ratio H/r 
which is zero at the disc inner edge and which increases as the magnetic support decreases. The disc is squeezed by 
the magnetic field if the magnetic pressure dominates over the thermal pressure. If the sound speed is about a tenth 
of the Keplerian velocity, then the magnetic and thermal pressures become comparable when the magnetic support 
is about 0.1. Therefore, we choose a profile of H/r which reaches a constant value (that we take to be 0.1) at the 
location where f3 = 0.1. 

We note that the vertical displacement of the disc at a location (r, ip) is: 

1l{Z z (r)} cos <p-l[Ur)] sin <p. (33) 



Figures [lj |3j show the magnetic support /3, the disc aspect ratio H/r and the surface density E/Ej versus w. Also 
displayed are Tl(^ z )/(26R ), i.e. £ Z /(26R ) for ip = versus x/R D , and -1(£ Z )/(28R ), i.e. £ >Z /{28R ) for tp = tt/2 
versus y/R , for a — 10~ 3 and 0.1. These quantities represent the disc vertical displacement within a factor 2SR a 
in the (ip — 0) and [<p — tt/2) half planes, respectively, or, equivalently, in the (x > 0, z) and (y > 0, z) half planes, 
respectively (see eq. |33|). The different figures correspond to different models. 

Figure [l] and ^ correspond to Aly's solution and wb = 0.0975. In Figure [l], 1//3 oc exp [100(n7 — 0.14)] + 1. The 
results do not depend significantly on the details of the H/r profile in the disc inner parts (we considered both the case 
where H/r increases almost linearly and exponentially with radius). We also checked that £ z hardly changes when tub 
varies between 0.09 and 0.099, providing we keep the same magnetic support (i.e. we change E accordingly). Figure]^ 
corresponds to (3 oc / {exp [40(-n7 — 0.14)] + 1}. The results are qualitatively the same as in Figure g. We also ran 
the case (3 oc l/ro 3 , which gave similar results, except that the positive values reached by £ z in that case were almost 
as large as the negative values. 

Figure ^ corresponds to Low's solution and wb = 0.09999999. Here again, the results hardly depend on the details 
of the profile of H/r in the disc inner parts. From equation ([27]), we see that £ 2 vanishes at w = vjb if E is non zero 
there. When wb is very close to Wi, we indeed check that £ z is almost zero at the disc inner edge. However, since B^ 
varies very rapidly near wb, Cz takes some finite value at the inner edge as wb is moved a little bit away from Wi. 
This is illustrated in Figure ||, where we plot f3 and £ 2 for different values of vjb, ranging from 0.09 to 0.09999999 (all 
the curves have the same parameters except for wb)- 

The calculations for the low viscosity a = 10 -3 illustrate the excitation of bending waves. Several oscillations 
corresponding to bending waves excited at the OLR are visible. These are damped by viscosity before they propagate 
very far. This damping is enhanced by the increasing wavenumber produced by the decreasing magnetic field. 
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For larger and probably more realistic viscosity a = 0.1, these waves are so heavily damped they are barely visible. 
The dominant displacement is near the disc inner edge in this case. 

We note that, when viscous damping is present, the response of the disc is not in phase with the perturbation, 
and the warp lags behind the dipole. That is, the disc is also twisted, and the twist is trailing. Mathematically, this is 
illustrated by the fact that the imaginary part of £ z is non zero, as seen in the different figures. It may seem surprising 
that, even when a is as low as 10~ 3 , the imaginary part of £ z is comparable in magnitude to its real part. This is 
because £ z varies rapidly with radius when a is small, and the imaginary part of £ z becomes important when |<9 2 £ z /3r 2 | 
is large (see eq. ]27]]). We checked however that the imaginary part of £ z becomes very small compared to its real part 
when a is reduced further. Practically, the existence of a twist means that the point of maximum elevation in the disc 
is not in the plane which contains the dipole axis and the rotation axis of the star, i.e. the (tp = 0), or equivalcntly 
(x, z), plane, but in a plane corresponding to a smaller (negative) tp. Or, in other words, the elevation does not vanish 
in the (tp — 7r/2), or equivalcntly (y, z), plane, but in a plane corresponding to a smaller value of tp. To illustrate this, 
we show in Fig. ^| the line of nodes (which joins the points of zero elevation) in the inner parts of the (x, y) plane for 
both Aly's and Low's models and a = 0.1. We clearly see that the line of nodes is trailing. If there were no twist, 
it would indeed coincide with the y-axis. The curves corresponding to different values of vjb or different models are 
very similar to those displayed in Fig. ||. For a = 10 -3 , the line of nodes is more tightly wrapped, because £ z oscillates 
more rapidly. 

In Fig. ^ we present a surface plot of the warped disc corresponding to the Low's model displayed in Fig. ^|. This 
clearly shows the presence of an obscuring wall like feature near the disc inner boundary and the trailing twist. 



5. Discussion and application to AA Tau 

The light curve of the classical T Tauri star AA Tau displays photometric, spectroscopic and polarimetric variations 



on timescales from a few hours to several weeks (Bouvier et al. 1999). The most striking feature of this light curve is 



a photometric variability with a period of 8 to 9 days, comparable to the expected rotation period of the star. This 



has been interpreted by Bouvier et al. (|1999|) as being due to the occultation of the star by a warp of the inner disc 



(the system is observed almost edge-on). They proposed that this warp be produced by the interaction of the disc 



with a stellar magnetic dipole tilted with respect to the disc rotation axis. Following Mahdavi & Kenyon (1998), they 
assumed that material at the disc inner edge, neglecting warping, would be more likely to accrete along the shorter 
field line than along the longer one connecting star and disc. In their model, based on this accretion geometry, the 
disc is elevated above the original equatorial plane at the location where the dipole axis is bent toward the disc. This 
geometry enables at least part of the stellar hot spots, located at the field line footpoints , to still be on the line of 
sight when occultation of the star is maximum, in agreement with the observations. In this model, the warp is not 
calculated self-consistently, and is not actually a real dynamical bending of the disc. 



The assumption by Mahdavi & Kenyon ( 1998| ) can become incorrect, precisely because it neglects the warping of 



the disc. If the inner disc is no longer in the original equatorial plane, then material at the inner edge no longer 'sees' 
the short and long field lines it would see if it were in the unperturbed disc plane. 



At the disc inner edge, since £l(Ri) = CIk(R%) — we see from equation (27) that the sign of £ z is the same as 
that of P' m (i.e., from eq. |l8[, the same as that of —B+B' r + ) for p = 0, and that determines the direction in which 
the disc bends. From Figure pi we see that for Low's solution the disc inner edge on the a;-axis is pushed below the 



equatorial plane, which is the opposite of what was anticipated by Mahdavi & Kenyon (1995). For Aly's solution, the 
disc is bent above the equatorial plane (see Fig. |l| and 0) . It is not surprising that the two cases give different results 
since changes sign from one case to another. 

To get the elevation of the disc above the equatorial plane, the real and imaginary parts of £ z have to be combined 
through equation (^). In Fig. 0, we show the maximum elevation t; Zt max (where the tp dependence has been taken into 
account in £ z ), normalized to unity, above the equatorial plane along the —p direction, versus — tp. This represents the 
maximum elevation of the disc material located in between the star and an observer looking at the disc almost edge-on 
and from above along the — tp direction. The curves correspond to Aly's and Low's models displayed in Fig. |l] and |[ 
respectively, and a = 0.1. Let us suppose that, to begin with, we look at the system along the tp — direction. Then, 
as the star rotates, our line of sight moves to smaller, negative values of tp (to the right along the x-axis of Fig. ^). 
In the case of Aly's model, this corresponds to the star being occulted at first and until it has rotated by about 90 
degrees. Then, as the star rotates further, it becomes less and less obscured, and it can be seen by the observer until 
it has completed a full rotation. In the case of Low's model, the star is not obscured at first. Occultation begins only 
after the star has rotated by about 180 degrees and lasts also for about a quarter of a period. We note that, depending 
on the direction of the line of sight, the radius at which the elevation is maximum varies (typically between 0.1 R 
and 0.2 R a ). It is unlikely that the observer would be able to know at which radius is the part of the disc responsible 
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for the occultation, as the whole structure rotates with the same frequency (that of the star). Indeed, it is important 
to stress again that a structure located at some radius r does not rotate with the local angular frequency but with 
that of the star (in other words, the warp appears steady in a frame rotating with the star). 

With Low's solution, the disc is pushed below the equatorial plane at the location where the dipole is bent toward 
the disc. Therefore the stellar hot spots, at the field line footpoints cannot be seen under condition of maximum 
occultation, unless they are relatively large. In Aly's case, at least part of the hot spots can still be observed at an 
occultation maximum because of the phase-lag between the disc response and the perturbation. Indeed, we see in 
Fig. ^ that, for the particular Aly's model represented there, the occultation is maximum when ip ~ —15°, and not 
when (p — 0. Therefore, it seems the observations favour Aly's model. 

To fit the observations, Bouvier et al. ( |1999 ) needed the amplitude of the warp to be about 0.3 times the disc inner 
radius. We see from Figures |j]-|3| that this is easily attained here. There, |£ z | / (2SRi) ranges from about 0.1 to 1. For S 
up to about 30° (for which shxS ~ 5 in radians), this gives |£ z | /Ri from 0.1 to 1. The lower value could even be made 
larger by considering a profile of H/r with an exponential rather than almost linear increase in Figure ||. Therefore 
only a moderate misalignment angle is required to produce |£ z | /Ri ~ 0.3. 

Figures show that when the viscosity in the disc is small (a = 10~ 3 for instance), the vertical displacement of 
the disc varies rapidly, on a scale smaller than the disc semi-thickness H. In that case, we expect the warp to become 
dispersive (Papaloizou & Lin 1995), and probably the disc to disrupt. It is likely that different situations arise when 
a disc interacts with a nonaligned dipole, depending on the detailed physics. In some cases we may get a moderate 
smoothly varying warp, as described above, in other cases it may be that the disc breaks up and reforms. This may 
also produce light curve variability. 

We note that the bending waves which are excited by the tilted dipole transport angular momentum (Papaloizou 
& Terquem 1995; Terquem 1998), which is deposited in the disc if the waves are damped. Since the dipole rotates 
faster than the disc in which the waves propagate, the resulting torque opposes accretion onto the star. This effect 
may produce additional variability, and will be the subject of another paper. 
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Fig. 1. Results for Aly's magnetic field. Top left panel: Magnetic support (3 and aspect ratio H/r for an almost 
linear {solid line) and exponential [dotted line) profile versus w = r/R a . Top right panel: Surface density S/S^ where 
Ej = versus zu = r/R Q Middle panels: lZ(t; z )/(2SR ), i.e. ^ Z /(28R ) for ip = versus x/R Q . Bottom panels: 

— T(£,z)/ (2<5 R ), i.e. £ Z /(25R ) for ip = ir/2 versus y/R - The different curves are for a = 10~ 3 (left panels) and 0.1 
(right panels ) and correspond to the almost linear (solid lines) and exponential (dotted lines) profile of H/r. Here 
zdb = 0.0975. Only the inner parts of the disc where the displacement is non zero are shown. Several oscillations 
corresponding to bending waves excited at the OLR are visible for a = 10~ 3 . These are damped by viscosity before 
they propagate very far. 
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Fig. 2. Same as Figure but for a different magnetic support (5 and surface density E. Only one profile of H/r is 
considered here. Results are similar, although the maximum of \£ z \ is now closer to the disc inner edge. 
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Fig. 3. Same as Figure but for Low's magnetic field with -cub = 0.09999999 and different magnetic support (3 and 
surface density S. The vertical displacement is almost constant and maximum very close to the disc inner edge, where 
it actually vanishes. 
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Fig. 4. Same as Figure |3|, but for different values of tub- The different curves correspond to tub = 0.09999999 (solid 
lines), 0.0999 (dotted lines), 0.099 (short-dashed lines), 0.0975 (long-dashed lines) and 0.09 (dotted-dahed lines). When 
tub is moved a little bit further away from the disc inner edge, the vertical displacement takes a finite value there. 
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Fig. 5. Line of nodes in the inner parts of the (x,y) plane for the same models as Figure ^ (Aly's magnetic field, filled 
triangles) and Figure |^ (Low's magnetic field, open triangles). Here the profile of H/r is almost linear and a = 0.1. 
We see that the line of nodes is trailing. If there were no twist, it would coincide with the y-axis. For a = 10~ 3 , it is 
more tightly wrapped. 
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Fig. 6. 3D view of the disc corresponding to the model represented in Figure y (Low's magnetic field). The dipole, its 
axis and the cc-axis are also represented. The vertical scale has been amplified for clarity. 
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Fig. 7. Maximum elevation ^ z . m ax (where the ip dependence has been taken into account in £ 2 ), normalized to unity, 
above the equatorial plane along the —ip direction, versus —ip in degrees. This represents the maximum elevation of the 
disc material located in between the star and an observer looking at the disc almost edge-on and from above along the 
—ip direction. The curves correspond to the models shown in Figure ^ (Aly's magnetic field, solid line) and Figure |^ 
(Low's magnetic field, dotted line). Here the profile of H/r is almost linear and a = 0.1. Depending on the direction 
of the line of sight, the radius at which the elevation is maximum varies (typically between 0.1 R a and 0.2 R ). For 
Aly's model, maximum occulation occurs for if in the range [— 7r/2,0], whereas for Low's model it occurs for ip in the 
range [— 37r/2,— w] or, equivalently, [w/2,tv]. 



